
;;======================================================================
;; Script to create an IDL structure akin to that created by
;;   create_rb3_struct.pro, except filled with data from the BGPS.
;;   This structure is designed to be a drop-in replacement for use
;;   with the IRDC foreground fractions modeling.
;;
;; NOTE: This script requires the UBC GLIMPSE data products to be
;;   located in the ./local/ubc_glimpse_proc/ directory.
;; ALSO: Requires the BGPS postage-stamp images located in the
;; ./local/postage/ directory.
;;
;;   **Run from the distance-omnibus directory.
;;     To run: create_bgps_rb3_struct [,isuff][,/DORAND]
;;======================================================================
;;   Changed to PRO on 07/04/12 to ease use.
;;   09/05/12, TPEB, Added IRAC aperture correction factor of 0.737
;;   for extended emission in the BAND 4 images.
;;   09/06/12, TPEB, Corrected the scattered light term introduced
;;   yesterday.
;;   10/19/12, TPEB, Moved IRAC correction term to
;;   process_glimpse_irdc.pro to make postage stamp images as close to
;;   correct as possible.


PRO CREATE_BGPS_RB3_STRUCT, isuff, DORAND=dorand
  
  IF ~ exist(isuff) THEN isuff = ''
  
  ;; Read in the BGPS distance database
  s = read_bgps_csv('bgps_distance_database.csv',csv,/ver)
  ppbeam = 23.8027594254d       ; # of BGPS pixels per beam
  
  ;; Read in Cara's temperature maps & get useful information
  TMAP_DIR = 'irdc_dist_model/battersby/'
  tmap = readfits(TMAP_DIR+'gaussfit_iter_beta175_temp_l30_conv25.fits',$
                  thdr,/silent)
  tlab = readfits(TMAP_DIR+'gaussfit_iter_beta175_column_l30_'+$
                  'conv25_label.fits',lhdr,/silent)
  terr = readfits(TMAP_DIR+'gaussfit_iter_beta175_temperr_l30_conv25.fits',$
                  ehdr,/silent)
  plot_wcs_axes,thdr,10,lonarr=larr
  lrange = minmax(larr)
  extast,thdr,tastr
  tsz = size(tmap, /DIM)
  
  ;; Next, read in the BGPS velocity structure, to extract VLSR
  restore,'local/bgps_velocity_struct.sav',/ver
  
  ;; Figure out how many objects have a dense gas composite spectrum
  dgind = WHERE( total(velocity.densegas,1) NE 0, ndg)
  velocity = velocity[dgind]
  
  ;; Check the UBC_PROC files for a list of filenames
  UBC_DIR = './local/ubc_glimpse_proc/'
  ufn = FILE_SEARCH(UBC_DIR+'bgps????.fits',COUNT=n_ubc)
  ubc_cnum = long(strmid(ufn,STRPOS(ufn[0],'bgps')+4,4))
  
  ;; Get the overlap between dense gas tracers and the UBC GLIMPSE proc
  ind = WHERE_ARRAY(velocity.cnum, ubc_cnum, n)
  ufn = ufn[ind]
  
  ;; Work backwards...
  ind2 = WHERE_ARRAY(ubc_cnum, velocity.cnum, n2)
  velocity = velocity[ind2]
  IF n NE n2 THEN $
     message,"Take a shot of gin, then figure out the size mismatch."
  
  print,'Total sample is '+string(n,format="(I0)")+' objects.'
  
  ;; Create RB3 structure for the BGPS sources
  rb3 = replicate(create_struct('name','','cl',0.,'cb',0.,'core','','cnum',0L,$
                                'dist',0.,'l',0.,'b',0.,'vlsr',0.,'s_peak',0.,$
                                'I_mir',0.,'I_min',0.,'del_Imir',0.,$
                                'ir_dx',0.,'ir_dy',0.,'ir_l',0.,$
                                'ir_b',0.,'ir_pa',0.,'C_meas',-1000.,$
                                'batt_td',0.,'batt_te',0.,$
                                'vlsr_flag',0b,'s_err',0.,$
                                'sig_imin',0.,'sig_imir',0.,'sig_C',0.,$
                                'aimir',0.,'aimin',0.,'asmir',0.,'asmin',0.,$
                                'ac',0.,'asc',0.,'I_min2',0.,'sig_imin2',0.,$
                                'grs',0b,'grs_match',0b,'grs_clemdist',0.,$
                                'grs_dist',0.,'grs_l',0.,'grs_b',0.,'grs_v',0.,$
                                'fitdist',0.d,'morph_d',0.d,'morph_width',0.,$
                                'morph_fwhm',0.,'morph_np',0.,'equal_prob',0b,$
                                'morph_sigm',0.,'morph_sigp',0.,'morph_pr',0.,$
                                'constrain',0b,'morph_constrain',0b), n)
  
  ;;====================================================================
  ;; Now, loop through each potential source...
  FOR i=0L, n-1 DO BEGIN
     
     IF ((i+1) MOD 100) EQ 0 THEN message,'Doing object #'+$
                                          string(i+1,format="(I4)")+' of '+$
                                          string(n,format="(I4)"),/inf
     
     ;; Since velocity and ubc are now the same size, loop through
     ;; velocity, since its cnum are in proper order
     j = where(ubc_cnum EQ velocity[i].cnum, nj)
     k = where(s.cnum EQ velocity[i].cnum, nk)
     IF nj EQ 0 THEN message,$
        "Take a shot of rum, then figure out why this object isn't in ubc."
     
     ;; Start filling in the rb3 structure
     rb3[i].name = s[k].name
     rb3[i].core = string(s[k].cnum,format="(I4)")
     rb3[i].cnum = s[k].cnum
     rb3[i].l    = s[k].glon_peak
     rb3[i].b    = s[k].glat_peak
     
     ;; Check on the longitude range -- GLIMPSE overlap!
     IF (rb3[i].l GT 65.25) THEN GOTO,ENDOFLOOP
     ;; IF (rb3[i].l LT 21.) THEN GOTO,ENDOFLOOP
     
     ;; Check latitude -- needs to be within safe boundaries of GLIMPSE
     IF ( abs(rb3[i].b) GT 1.) THEN GOTO,ENDOFLOOP
     
     ;; INSERT PROPER FLUX MEASUREMENT HERE -- structure expects mJy
     ;; Bolocat fluxes x 1.5+-0.15
     rb3[i].s_peak = s[k].flux_40 * 1.d3 * 1.5d
     sig_ss = (s[k].err_flux_40 > 0.015) * 1.d3 * 1.5d
     rb3[i].s_err  = sqrt(sig_ss*sig_ss + rb3[i].s_peak*rb3[i].s_peak*0.0225d)
     
     ;; Update VLSR with velocity.densegas, if available, otherwise use
     ;; GRS 13CO
     IF total(velocity[i].densegas) NE 0 THEN BEGIN
        yfit = GAUSSFIT(velocity[i].v_std,velocity[i].densegas, A, NTERMS=3)
        rb3[i].vlsr = A[1]
        rb3[i].vlsr_flag = 1b
     ENDIF ELSE BEGIN
        yfit = GAUSSFIT(velocity[i].v_std,velocity[i].grs13co, A, NTERMS=3)
        rb3[i].vlsr = A[1]
        rb3[i].vlsr_flag = 2b
     ENDELSE
     
     ;; Check on the longitude range -- can't use objects w/in 21
     ;;                                 deg of GC UNLESS part of
     ;;                                 Molecular Ring structure
     IF (rb3[i].l LT 21.) THEN $
        IF (rb3[i].vlsr GT (40./15.3*(rb3[i].l-7.5)+40.)) || $
        (rb3[i].vlsr LT (30./15.3*(rb3[i].l-7.5))) THEN GOTO,ENDOFLOOP

     ;; Everything else from the structure!  We are using GLIMPSE data
     ;; products produced by the script irdarkcloudsfinal3.pro, modified
     ;; from the one produced by the UBC group under Erik. 
     
     ;; String of cnum
     scnum = string(s[k].cnum,format="(I04)")
     imir_img = readfits(UBC_DIR+'Imir'+scnum+isuff+'.fits',$
                         imir_hdr,/silent)
     nons_img = readfits(UBC_DIR+'nonsmooth'+scnum+isuff+'.fits',$
                         nons_hdr,/silent)
     smoo_img = readfits(UBC_DIR+'smooth30arc'+scnum+isuff+'.fits',$
                         smoo_hdr,/silent)
     
     ;; For postage-stamp images going off the edge of the GLIMPSE
     ;; map, replace 0's with NANs -- which can more easily be
     ;; dealt with later.
     inan = WHERE(nons_img EQ 0, nnan)
     IF nnan GT 0 THEN BEGIN
        imir_img[inan] = !values.f_nan
        nons_img[inan] = !values.f_nan
        smoo_img[inan] = !values.f_nan
     ENDIF
     
     ;; TO CALCULATE I_MIR & I_min (with their uncertainties), we use
     ;;    the 40" Bolocat aperture (to align with BGPS flux
     ;;    measurement).
     ;; Extract pixel information using the 40" aperture method
     ;;    from Bolocat
     wtmask  = aperture_wtmask(smoo_hdr, 40., rb3[i].l, rb3[i].b)
     ind     = where(wtmask GT 0.5)
     wtmask2 = aperture_wtmask(imir_hdr, 160., rb3[i].l, rb3[i].b)
     ind2    = where(wtmask2 GT 0)
     
     
     ;; Calculate the extracted number
     ;; rb3[i].I_min2 = total(wtmask*smoo_img, /NAN) / total(wtmask)
     ;; I_min should be the minimum pixel value within this aperture,
     ;; since we have already smoothed it.
     rb3[i].I_min = min( smoo_img[ind], /NAN )
     ;; Actually, I_min should be the aperture weighted average of the
     ;; non-smoothed map, since this sidesteps the dumping of bright
     ;; into dark areas.
     rb3[i].I_min2 = total(wtmask*nons_img, /NAN) / total(wtmask)
     
     
     rb3[i].I_MIR = total(wtmask2*imir_img,/NAN) / total(wtmask2)
     
     ;; Calculate the uncertainty
     V1 = total(wtmask,/NAN)
     V2 = total(wtmask*wtmask,/NAN)
     rb3[i].sig_imin2 = sqrt( total(wtmask*(smoo_img - rb3[i].I_min2)^2,$
                                    /NAN) * (V1 / (V1*V1 - V2)) )
     rb3[i].sig_imin = stddev( smoo_img[ind], /NAN )
     
     
     V1 = total(wtmask2,/NAN)
     V2 = total(wtmask2*wtmask2,/NAN)
     rb3[i].sig_imir = sqrt(total(wtmask2[ind2]*$
                                  (imir_img[ind2]-rb3[i].I_mir)^2,/NAN) * $
                            (V1 / (V1*V1 - V2)) )
     
     ;; Calculate contrast
     rb3[i].C_meas = ( 1.d - (rb3[i].I_min / rb3[i].I_mir) )
     rb3[i].sig_C  = sqrt( (rb3[i].i_min*rb3[i].i_min/$
                            rb3[i].i_mir/rb3[i].i_mir)*$
                           (rb3[i].sig_imir*rb3[i].sig_imir/$
                            rb3[i].i_mir/rb3[i].i_mir + $
                            rb3[i].sig_imin*rb3[i].sig_imin/$
                            rb3[i].i_min/rb3[i].i_min) )
     
     ;; Calculate ALTERNATE (i.e. original) values for GLIMPSE
     rb3[i].aimin = rb3[i].i_min
     rb3[i].aimir = MEAN( imir_img, /NAN )
     rb3[i].asmin = 0.3
     rb3[i].asmir = ( MAX(imir_img,/NAN) - MIN(imir_img,/NAN) ) / $
                    MIN(imir_img,/NAN) * rb3[i].aimir
     rb3[i].ac  = ( 1.d - (rb3[i].aimin / rb3[i].aimir) )
     rb3[i].asc = sqrt( (rb3[i].aimin*rb3[i].aimin/rb3[i].aimir/rb3[i].aimir)*$
                        (rb3[i].asmir*rb3[i].asmir/$
                         rb3[i].aimir/rb3[i].aimir + $
                         rb3[i].asmin*rb3[i].asmin/$
                         rb3[i].aimin/rb3[i].aimin) )
     
     
     ;; Check for correspondence with Cara's temperature maps from Hi-GAL
     AD2XY, rb3[i].l, rb3[i].b, tastr, x, y
     IF (x GE 0) && (x LT tsz[0]) && (y GE 0) && (y LT tsz[1]) THEN BEGIN
        
        ;; Next, check to see if location is within studied region
        IF tlab[x,y] THEN BEGIN
           
           rb3[i].batt_td = tmap[x,y]
           rb3[i].batt_te = terr[x,y]
           
           str = rb3[i].name+'_'+rb3[i].core+'  Td = '+$
                 string(tmap[x,y],terr[x,y],FORMAT="(F4.1,' +- ',F4.1)")+$
                 ' K   LABEL = '+$
                 string(tlab[x,y],FORMAT="(I0)")
           print,str
           
        ENDIF
     ENDIF
     
     ENDOFLOOP:
  ENDFOR
  
  ;; Save full list for comparisons, as needed
  save,rb3,filename='./irdc_dist_model/bgps_rb3_uncut.sav',/ver
  
  
  ;; Check for POSITIVE IRDC contrast (physical concordance with the model)
  ;; Check for sources with suspicious VLSR ( < -30 km/s)
  goodind = WHERE( (rb3.c_meas GT 0.01) AND (rb3.vlsr GE -30.), ngood)
  rb3 = rb3[goodind]
  
  print,'Number of useable objects: ',n_elements(rb3)
  
  ;; Check against BAD STARSUB list
  readcol,'./irdc_dist_model/bad_starsub_sources.txt',bad_ss_cnum,format='I'
  good_ss_ind = MISSING(bad_ss_cnum, rb3.cnum, n_good_ss)
  print,'# Good: ',n_good_ss
  rb3 = rb3[good_ss_ind]
  
  ;;====================================================================
  ;; We now have the final list of sources!
  ;; Save the structure
  save,rb3,filename='./irdc_dist_model/bgps_rb3'+isuff+'.sav',/ver
  
  spawn,'cp ./irdc_dist_model/bgps_rb3'+isuff+$
        '.sav /d4/idl_lib/distance-omnibus/irdc_dist_model/.'
  
  ;;=============================================================
  ;; For BASIC MODEL ONLY...
  IF isuff EQ '' THEN BEGIN
     
     ;;==================================
     ;; Sources with temperatures from Cara's work...
     ind = where(rb3.batt_td NE 0., ntd)
     print,ntd
     rb3=rb3[ind]
     save,rb3,filename='./irdc_dist_model/bgps_rb3_temp.sav',/ver
     
     spawn,'cp ./irdc_dist_model/bgps_rb3_temp.sav '+$
           '/d4/idl_lib/distance-omnibus/irdc_dist_model/.'
     
     ;; Call the testing set creation routine
     CREATE_TESTING_RB3_SETS, DORAND=dorand
     
  ENDIF
  
  ;; Clear out the memory a liitle
  undefine,velocity,s,ubc,rb3
  
END
